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ABSTRACT 

The stability of radially anisotropic spherical stellar systems in modified Newtonian 
dynamics (MOND) is explored by means of numerical simulations performed with the 
iV-body code n-mody. We find that Osipkov-Merritt MOND models require for stabil- 
ity larger minimum anisotropy radius than equivalent Newtonian systems (ENSs) with 
dark matter, and also than purely baryonic Newtonian models with the same density 
profile. The maximum value for stability of the Fridman-Polyachenko-Shukhman pa- 
rameter in MOND models is lower than in ENSs, but higher than in Newtonian models 
with no dark matter. We conclude that MOND systems are substantially more prone 
to radial-orbit instability than ENSs with dark matter, while they are able to sup- 
port a larger amount of kinetic energy stored in radial orbits than purely baryonic 
Newtonian systems. An explanation of these results is attempted, and their relevance 
to the MOND interpretation of the observed kinematics of globular clusters, dwarf 
spheroidal and elliptical galaxies is briefly discussed. 
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1 INTRODUCTION 

Modified Newtonian dynamics (MOND) was originally pro- 
posed by Milgrom (1983) to explain the rotation curves of 
disk galaxies without invoking the presence of dark matter 
(DM) and, over the years, it has been successful at reproduc- 
ing the observed kinematics of several galaxies (e.g. Sanders 
& McGaugh 2002; Sanders & Noordermeer 2007; Swaters, 
Sanders, & McGaugh 2010, and references therein). How- 
ever, to some extent, the MOND and DM interpretations of 
the kinematics of galaxies can be degenerate. For instance, 
a MOND rotation curve can be also described in the Newto- 
nian gravity invoking a DM halo such that the total gravi- 
tational field is the same as the MOND one. More generally, 
given a MOND system, it is possible to construct the equiv- 
alent Newtonian system (ENS), i.e. the Newtonian system 
with DM in which the visible matter has the same phase- 
space distribution as in the MOND system (Milgrom 1986, 
2001; Nipoti et al. 2007b; Nipoti, Londrillo & Ciotti 2007c; 
Nipoti et al. 2008). It should be noted, however, that the 
physical viability of the ENS is not guaranteed, as for some 
configurations the density of the associate DM halo turns 
out to be negative (Milgrom 1986). Though a MOND sys- 
tem and its ENS are, by construction, indistinguishable from 
a kinematic point of view (for instance, as far as the rota- 
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tion curve^ or the velocity-dispersion profile are concerned), 
in general they are not identical from a dynamical point of 
view: for instance, two-body relaxation, dynamical friction 
and galaxy merging act differently in the two cases (Ciotti 
& Binney 2004; Nipoti et al. 2007c, 2008). 

As already recognized, MOND and ENSs might differ 
substantially also in terms of stability. As one of the original 
motivations for invoking DM halos in disk galaxies was that 
purely baryonic Newtonian disks are prone to bar-like insta- 
bility (Ostriker & Peebles 1973), it is not surprising that the 
study of dynamical stability in MOND has focused mainly 
on disks (Milgrom 1989; Christodoulou 1991; Brada & Mil- 
grom 1999). Here we consider instead the so-called radial- 
orbit instability, which is relevant to pressure-supported stel- 
lar systems. As well known, in the context of Newtonian 
gravity the amount of radial orbits in stellar systems is lim- 
ited not only by the requirement of phase-space consistency 
(i.e. positivity of the distribution function; see, e.g., Ciotti & 
Pellegrini 1992; An & Evans 2006; Ciotti & Morganti 2009, 
2010a,b) but also by the fact that very radially anisotropic 
spherical systems are unstable (Henon 1973; Polyachenko & 
Shukhman 1981; Fridman & Polyachenko 1984; Merritt & 

^ The MOND rotation curve of a disk galaxy can be always re- 
produced with a spherical DM halo and Newtonian gravity. How- 
ever, the vertical kinematics differs in the two cases (Nipoti et 
ah 2007b), because the DM halo of the ENS of a disk galaxy is 
non-spherical. 
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Aguilar 1985; Barnes, Hut, & Goodman 1986; May & Bin- 
ney 1986; Palmer & Papaloizou 1987; Dejonghe & Merritt 
1988; Saha 1991; Weinberg 1991; Berlin et al. 1994; Hjorth 
1994; Meza & Zamorano 1997; Trenti & Berlin 2006; Barnes, 
Lanzel, & Williams 2009). In Newtonian gravity DM halos 
tend to have a stabilizing effect against radial-orbit instabil- 
ity (Stiavelli & Sparke 1991; Meza & Zamorano 1998; Nipoti, 
Londrillo, & Ciotti 2002), which suggests that MOND sys- 
tems might be more prone to this kind of instability than 
their ENSs. If confirmed, this could provide a discriminant 
between the two theories when interpreting the velocity dis- 
persion profiles in pressure-supported stellar systems. For 
these reasons, in this work we use A^-body simulations to 
explore the stability of radially anisotropic MOND spheri- 
cal systems and, for comparison, of their ENSs and of purely 
baryonic Newtonian systems with the same density distri- 
butions and anisotropy profiles. 

The paper is organized as follows. In Section 2 the 
galaxy models used in the simulations are described. Their 
phase-space consistency is discussed in Section 2.1 and their 
stability in Section 3. Section 4 concludes. 



2 GALAXY MODELS 

We consider MOND as modified gravity in the non- 
relativistic formulation^ of Bekenstein & Milgrom (1984), 
in which the Poisson equation V^(^^ — AnGp is replaced by 



ao 



= AttGp, 



(1) 



where 



is the standard Euclidean norm, 

-10, 



is the MOND 

gravitational potential, ao — 1-2 x 10~^"ms~'^ is a char- 
acteristic acceleration, and jj, is the so-called interp olating 
function (in the present work we adopt jj.{y) — y/ yj'l + y^; 
Milgrom 1983). The MOND gravitational acceleration is g = 
—Vcf), just as the Newtonian acceleration is = — V^'^. For 
a system of finite mass, — >■ as ||x|| — >■ oo, where x is 
the position vector. 

As well known, from the Poisson equation and equa- 
tion (1) it follows that the MOND and Newtonian gravita- 
tional accelerations are related by n{g/ao) g = g'^-l-S, where 
g = ||g|| , and S is a solenoidal field dependent on the specific 
p considered. Since in general S 7^ 0, and its expression is un- 
known a prion, standard Poisson solvers cannot be used to 
develop MOND A-body codes, in which equation (1) must 
be solved at each time step (see Brada & Milgrom 1999; 
Ciotti, Londrillo & Nipoti 2006; Tiret & Combes 2007). In 
the present work we use our original MOND A-body code N- 
MODY (Nipoti, Londrillo & Ciotti 2007a; Londrillo & Nipoti 
2009, see Section 3.1). 

The initial conditions of the simulations are A-body re- 
alizations of galaxy models with the stellar component de- 
scribed by a spherical 7-model (Dehnen 1993; Tremaine et 
al. 1994) with density distribution 



p,{r) = 



■ 7 



4n r~<{ri, + r] 



4 — , ' 



(2) 



^ Recently an alternative non-relativistic formulation of MOND, 
dubbed QUMOND, has been proposed (Milgrom 2010). 



where M, is the total stellar mass, r* is the scale radius, 
and 7 is the negative of the inner logarithmic density slope 
(0 < 7 < 3). For simplicity we restrict to the cases 7 = 
and 7 = 1 (Hernquist 1990); we recall that for 7 7^ 2 the 
Newtonian potential is 



(r) = 



GM, 



r*(2-7) 



r*+r 



2-7 



1 



(3) 



For these models the MOND potential, required to dis- 
tribute the particles in the velocity space, is easily calcu- 
lated from the Newtonian one, as in spherical symmetry 
S = 0. The ENS associated with a model of stellar density 
p* and MOND potential </> has a total density V^0(r)/47rG. 
However, the resulting DM halo density (obtained after sub- 
traction of p*) in principle may be negative or have a non- 
monotonic radial trend. Fortunately, it can be proved that 
the DM halos of the ENS derived from the spherical 7 = 
and 7 = 1 models are everywhere positive. Instead, possible 
non-monotonicity of the DM density distribution makes the 
consistency of the halo a non-trivial request, as we will dis- 
cuss in Section 2.1 (see also Ciotti, Morganti, & de Zeeuw 
2009). 

In order to impose a tunable amount of radial 
anisotropy on the initial conditions, we adopt the widely 
used Osipkov-Merritt (OM) parametrization (Osipkov 1979; 
Merritt 1985), in which the stellar distribution function de- 
pends on energy and angular momentum per unit mass 
through the variable Q = £ — J^/2rf , where £ = ^ — v'^ /2 is 
the relative energy, « = ||v|| is the velocity modulus, ^ — —(j> 
is the relative potential, J is the angular momentum modu- 
lus per unit mass, and ra is the so-called anisotropy radius. 
The anisotropy parameter (e.g. Binney & Tremaine 2008) 
is /3(r) = r^/(ra + r'^): for r 2> Ta the velocity dispersion 
tensor is radially anisotropic, while for r <^ the tensor 
is nearly isotropic. In the limit ra — >■ 00, Q = f and the 
velocity dispersion tensor becomes globally isotropic. 

In the purely baryonic Newtonian models, the distribu- 
tion function can be written as 
cQ 



/n(Q) 



1 



VSTV^dQj^ d*?^/QT^' 



(4) 



where = —(p^ is the Newtonian relative potential and 



g,{r) = l + _ p,(r) 



(5) 



is the so-called augmented density (see, e.g., Binney & 
Tremaine 2008). 

In the MOND (and in the ENS) cases a similar formula 
holds, i.e. 

rQ 



/m(Q) = 



1 d 



dg, d* 



VStv^ dQ d* 70^' 



(6) 



where now ^ is the MOND relative potential that, by con- 
struction, is also the total potential of the ENS. Note that, at 
variance with equation (4) , the lower integration limit is now 
—00: due to the far-field logarithmic behaviour of the MOND 
potential in isolated systems of finite mass, the distribution 
function must be positive for all values of Q < in analogy 
with finite-mass Newtonian system in a total isothermal po- 
tential, where stars of all energies are bound (e.g. Ciotti et 
al. 2009). Apart from a few exceptions (see Section 3.2.3), in 
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Table 1. Main properties of the studied families of models. 



Family 


Gravity 


(AfDM/A4'.)half 


»'ac/r-half 


'•as/fhalf 




Shalf.s 


(1) 


(2) 


(3) 


(4) 


(5) 


(6) 


(7) 


No 


Newton 





0.092 


~0.8 


~1.75 


~1.34 


MqKio 


MOND 





0.092 


~1.0 


~2.33 


~1.26 


EqKio 


Newton 


0.76 


0.092 


~0.5 


~3.48 


~1.83 


Mqko.oi 


MOND 





0.092 


~1.0 


~2.48 


~1.28 


Eqko.oi 


Newton 


50.1 


0.092 


<0.2 


>7.48 


>4.06 


Ni 


Newton 





0.053 


-0.8 


~1.64 


~1.28 


MiKio 


MOND 





0.053 


~0.8 


~2.32 


~1.31 


EiKio 


Newton 


0.31 


0.053 


~0.6 


~2.70 


~1.50 


MiKo.oi 


MOND 





0.053 


~0.9 


~2.61 


~1.31 


EiKo.oi 


Newton 


32.7 


0.053 


<0.2 


>6.86 


>3.64 



(1) Name of the family of models: N are the Newtonian purely baryonic models, M the MOND models, and E the ENS models. The 
first subscript is the value of the parameter 7 identifying the stellar density profile (equation 2); the second subscript is the value of the 
internal acceleration ratio k = GMt /aorl. (2) Gravity law. (3) DM-to-stellar mass ratio within the stellar half-mass radius rj^aif 
(riiaif — 3.84r« for 7 = and rji^if — 2.41r« for 7 = 1). (4) Normalized critical anisotropy radius for consistency, as determined from 
the GDSAI. Models with ra < rac are necessarily inconsistent. (5) Normalized minimum anisotropy radius for stability. (6) Maximum 
value for stability of the Fridman-Polyachenko-Shukhman parameter. (7) Maximum value for stability of the radial-to-tangential kinetic 
energy ratio measured within rj^^if. In the columns (5-7), upper limits on ra and lower limits on and ^half a mean that no instability 
has been detected in the most anisotropic case explored for that family. 



the simulations of the ENSs the halo is maintained "frozen" 
(i.e. it acts as a fixed external potential), so its distribution 
function is not needed. 

2.1 Consistency 

As well known, in (single- or multi-component) OM models 
it is possible to determine a critical value of the anisotropy 
radius of each density component, Tac, so that for < r^c 
the models are inconsistent, i.e., f{Q) < for some admis- 
sible value of Q. A necessary condition for consistency of 
OM models is dg,{r)/dr < at all radii (Ciotti & Pelle- 
grini 1992; Ciotti 1996, 1999; Ciotti et al. 2009). As recently 
proved in Ciotti & Morganti (2010a, b; see also Van Hese, 
Baes &i Dejonghe 2011; An 2011) this condition is rigorously 
equivalent to the Global Density Slope- Anisotropy Inequal- 
ity (GDSAI), i.e. to the requirement that at each radius the 
negative of the logarithmic slope of the stellar density dis- 
tribution cannot be less than twice the local value of the 
anisotropy parameter, 7*(r-) > 2/?(r-). Remarkably, it is easy 
to show that the GDSAI is independent of the gravity law 
holding the system together, so that it holds not only for 
Newtonian multi-component systems (such as the ENS), but 
also in the MOND case^. 

We applied the GDSAI to our families of models, deter- 
mining for each stellar density profile the value of rac (see 
Table 1), and the obtained limits coincide with those deter- 
mined in Ciotti (1999). As these are just necessary condi- 
tions for consistency, the positivity of the distribution func- 
tion must still be checked numerically for ra > rac- We also 
note that from the inequality dQt,{r)/dr < it follows that 
in case of isotropy (ra = cxj) a spherical density distribution 
must be necessarily monotonically decreasing for increasing 
r, and this imposes quite strong constraints on the DM halos 

^ This however is not true for the sufficient conditions for con- 
sistency (see Ciotti & Morganti 2010a, b) 



of physically admissible ENS. If the halo has a central "hole" 
(which is not unusual for ENSs; see Nipoti et al. 2007c), it 
can be physically realized only if its orbital distribution is 
tangentially biased, at least in the internal regions. In prac- 
tice, this is not an issue in simulations with frozen DM halos, 
but it is an important problem when setting up initial con- 
ditions for simulations in which the DM halo is "live" (i.e. 
modelled with particles; see Section 3.2.3). 



3 STABILITY 

We now describe the main results of the A'^-body simulations. 
Before, we briefly illustrate the numerical code n-mody and 
the tests performed. 

3.1 The numerical simulations 

We ran A'^-body simulations with n-mody, our original par- 
allel three-dimensional particle-mesh code that can be used 
to follow the evolution of either MOND or Newtonian col- 
lisionless stellar systems (Nipoti et al. 2007a; Londrillo & 
Nipoti 2009). In previous papers we have already used N- 
MODY to demonstrate signiflcant differences in violent relax- 
ation, merging and dynamical friction in MOND and New- 
tonian dynamics (Nipoti et al. 2007a,c; Ciotti et al. 2007; 
Nipoti et al. 2008) . We refer the readers to these papers and 
to Londrillo & Nipoti (2009) for a more detailed description 
of the code. In the present study the spherical grid has 64 ra- 
dial nodes, 32 nodes in colatitude and 64 nodes in azimuth 
ip, and the total number of particles is A'part — 8 x 10^. We 
verified with convergence experiments that these numbers 
of particles and grid points exclude important discreteness 
effects. In particular we reran some of the simulations with 
a grid 128 x 64 x 128 and A^part 6.4 x 10^, finding that 
the results are almost identical to the corresponding lower- 
resolution cases. 

In each simulation the initial conditions consist of an 
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Figure 1. Stellar (dotted curves) and DM (solid curves) density profiles of ENSs corresponding to spherical MOND 7 = (left column) 
and 7 = 1 (right column) models, for two different values of the dimensionless acceleration ratio k = GMt /cio''*- 



A''-body realization of the stellar distribution of an isolated 
equilibrium galaxy model, in which the particles are dis- 
tributed in phase space with the standard rejection tech- 
nique, using the numerically recovered distribution func- 
tions. In the ENS simulations, even if the DM halo is frozen, 
the systems are truncated exponentially at a radius rt, 
so that the potential well is finite. As a rule, we adopt 
rt = lOr,, but we ran some of the simulations also with 
rt — lOOr* and rt = lOOOr,, finding that the stability prop- 
erties of a model do not depend on the specific choice of the 
truncation radius. 

Following Nipoti et al. (2007c), we identify each MOND 
initial condition by fixing a value of the dimensionless 
internal acceleration parameter k = GMt/agr'^, so Af» 
and r* are not independent quantities. In physical units, 
introducing the quantity M*,io = M^/W^^Mq, r, ~ 
3Ak,~^^^ My^Qkpc, the time and velocity units are = 
^rllGM, ~ 29.7K-=*/^M,\/oMyr, and = r,/t^ ~ 
112K^/*M^^''iokms"\ The simulations are evolved up to 
lOOtdyn with a time step /St = O.Olfdyn, where fdyn is the 
characteristic dynamical time of the system. In the purely 



baryonic Newtonian 7-models, we adopt the standard half- 
mass radius value 

_ / Stt _ 7ri, 

= V 16Gphaif " V2 [21/(3-7) _ ' 

where phaif = 3M,/87rrjJjjif is the mean stellar density inside 
the half-mass radius rhaif (we recall that rhaif — 3.84r, for 
7 = and rhaif — 2.41r, for 7 = 1). In the corresponding 
MOND models (and their ENSs) , tdyn is given by the above 
expression multiplied by /v^, where is the circular 
velocity at rhaif for the purely baryonic Newtonian system 
and is the same quantity for the MOND system. 

Following Nipoti et al. (2002), in order to determine 
whether a given model is unstable, we check its departures 
from spherical symmetry by monitoring the evolution of its 
intrinsic axis ratios c/a and b/a (where a, b and c are the 
longest, intermediate and shortest axes of the inertia el- 
lipsoid of the stellar distribution). As preliminary tests we 
ran simulations in which the initial conditions are isotropic 
purely baryonic Newtonian, MOND and ENSs. We found 
that numerical uncertainties (due to the finite number of 
particles) on c/a and b/a never exceed 1% over lOOtdyn 
of evolution of these stable isotropic models. As a conse- 
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= MOND (M) O = Purely baryonic Newtonian (N) <^ = Equivalent Newtonian (E) 
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Figure 2. Final axis ratio (c/a)^^ as a function of the initial anisotropy radius normalized to the half-mass radius for MOND (stars), 
purely baryonic Newtonian (circles), and ENSs (squares) with different values of 7 and k. In each panel the horizontal dashed line marks 
the fiducial threshold value for stability, (c/a)^^ = 0.99; the arrows indicate the range of r-a/rjj^if corresponding to stable models. 



quence, we define unstable the models for which the axis 
ratio (c/a)fln after 100 t^yn is smaller than a fiducial thresh- 
old value c/a — 0.99. In the simulations the onset of the 
instability is just due to numerical noise produced by dis- 
creteness effects in the initial conditions. As expected, we 
found that an exact determination of the stability threshold 
for a given family of models is not straightforward: while for 
strongly anisotropic initial conditions the onset of the in- 
stability is apparent and the numerical models settle down 
into a final equilibrium configuration in a few dynamical 
times, for nearly stable initial conditions the instability can 
be characterized by very slow growth rates and its effects 
become evident even after tens of tdyn- 

3.2 Results 

For a given initial stellar density profile (i.e. 7 = or 7 = 
1 in equation 2), we explored, besides the purely baryonic 
Newtonian case (models No and Ni in Table 1), two families 
of MOND systems (and the corresponding families of ENSs), 
characterized by the values of the acceleration ratio n = 
0.01 and k = 10 (see Table 1 for a summary). The cases 



with K — 0.01 correspond to low-acceleration systems with 
internal accelerations everywhere much lower than ao, in the 
so-called deep-MOND regime. Their associated ENSs are 
therefore totally DM dominated systems (Table 1, Column 
3). On the other hand, the k = 10 cases are systems with 
relatively high accelerations within rhaif, corresponding in 
the Newtonian context to systems dominated by baryons in 
the central regions and by DM only at r>rhaif (Table 1, 
Column 3). The four panels of Fig. 1 show the baryonic and 
DM initial density profiles of the ENSs corresponding to the 
K = 0.01 and k = 10 MOND systems, for the two explored 
values of 7. As can be seen, the associated DM halos do not 
show a central hole and so, at least from this point of view, 
are not unrealistic (see Section 2.1). 

3.2.1 Minimum value for stability of the anisotropy radius 

For each family of models we have a set of eight simula- 
tions, in which the initial conditions differ only in the value 
of the normalized anisotropy radius ra/rhaif- This is appar- 
ent in Fig. 2, where we plot the final axis ratio (c/a)fln vs. 
the initial value of ra/rhaif for all the explored models of 
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Figure 3. Final axis ratio (c/a)(in as a function of the Fridman-Polyachenko-Shukhman parameter ^, for the same models as in Fig. 2. 
The arrows indicate the range of ^ corresponding to stable models. 



the ten families. The stability properties of the models can 
be inferred directly from Fig. 2: it is clear that the value of 
Ta at which instability appears is largest for MOND mod- 
els, smallest for the corresponding ENSs, and intermediate 
for purely baryonic Newtonian systems. Because of the very 
nature of the radial-orbit instability (very slow growth rates 
for marginally unstable systems) , a precise determination of 
the value of the minimum anisotropy radius for stability ras 
is difficult. What can be estimated robustly with A'^-body 
simulations is a fiducial value of ras separating apparently 
unstable systems from bona fide stable systems. Adopting 
as fiducial threshold for stability {c/a)fin = 0.99 (horizontal 
dashed line in each panel of Fig. 2), we estimate for each 
family of models the minimum anisotropy radius for stabil- 
ity ras, which is reported (in units of rhaif) in Table 1. 

As expected, the differences in the stability limit be- 
tween MOND and ENSs are more evident in systems with 
lower values of the acceleration ratio k (or, from a Newto- 
nian point of view, for more DM dominated systems). For 
instance, if we consider 7 = models with k = 0.01, which 
are representative of low-surface density systems with flat 
inner stellar density profile (such as dwarf spheroidal galax- 
ies), we find ras/rhait ~ 1 for MOND and ras/rhaif < 0.2 



for the ENS; very similar results are obtained for the more 
peaked Hernquist 7=1 models. Unsurprisingly, Newtonian 
models with the DM halo are less subject to the radial-orbit 
instability than the purely baryonic Newtonian models. 

Thus, from this first analysis we conclude that, when 
using ra/rhait as an indicator of the amount of admissible 
radial orbits, MOND systems are more prone to radial-orbit 
instability than their associated ENSs, and also than cor- 
responding purely baryonic Newtonian models. However, as 
we will see in the next Section, the comparison between the 
MOND and the associated purely baryonic Newtonian fami- 
lies is subtle: the fact that the latter typically admit smaller 
ra/rhaif values does not imply that (globally) they can sus- 
tain more radial kinetic energy. 

3.2.2 Maximum value for stability of the ^ parameter 

Stability limitations expressed in terms of ra are particu- 
larly relevant to observational works, as ra enters directly 
in the Jeans equations that are routinely used to fit the 
velocity-dispersion profiles of stellar systems. However, from 
a deeper point of view, the value of ra (loosely speaking, the 
radius outside which orbits are mainly radial) is not a ro- 
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bust measure of the fraction of kinetic energy that is stored 
in radial orbits, which, at least in the Newtonian context, 
is believed to be the main indicator of the tendency to de- 
velop radial-orbit instability. More specifically, it has been 
argued (Polyachenko & Shukhman 1981; Fridman & Poly- 
achenko 1984) that a proper criterion for stability for New- 
tonian self-gravitating systems can be expressed in terms of 
the global anisotropy parameter ^ = 2T-c/Tt, where Tr and 
Ti, = -\- Tip axe the radial and tangential components of 
the kinetic energy tensor, respectively. Global isotropy cor- 
responds to a value of the Fridman-Polyachenko-Shukhman 
parameter ^ = 1. Indications exist that there is a critical 
value such that only systems with ^ < are stable, and 
it is widely accepted that ~ 1.7 ± 0.25, relatively inde- 
pendent of the specific density distribution. 

A priori we do not have reasons to expect that ^ can 
be used as a stability indicator also in MOND. In any case, 
by definition measures how much radial anisotropy can 
be supported by a stable system, so it is interesting to dis- 
cuss the stability properties of our MOND, purely baryonic 
Newtonian and ENSs in terms of ^. The values of for the 
families of models studied in the present work are reported 
in Table 1, and in Fig. 3 we show all the models in the 
(c/a)fin-C plane to be compared with the analogous Fig. 2. 

For the purely baryonic Newtonian models the inter- 
pretation of these numbers is straightforward: we find ~ 
1.6 — 1.8, in agreement with the standard criterion and with 
previous numerical studies (see, e.g., Nipoti et al. 2002). For 
the MOND models we find 2.3<^s<2.6, with depending 
on both the stellar density profile and the acceleration ratio 
K. These values are substantially higher than those found 
for the corresponding purely baryonic Newtonian models, 
indicating that, for given density distribution, MOND sys- 
tems can sustain more radial kinetic energy than Newto- 
nian systems without DM. This is not in contrast with the 
findings reported in the previous Section, i.e. that MOND 
models have larger r^s than purely baryonic Newtonian mod- 
els. In fact, for fixed density profile and value, a MOND 
model has higher ^ than a purely baryonic Newtonian model, 
because more kinetic energy is stored in the outer parts, 
where orbits are radially biased and the gravitational field of 
MOND systems is stronger'*. Apparently, this effect compen- 
sates the larger values of admissible ra. Instead, the ENSs 
are again found able to sustain systematically higher values 
of than the corresponding families MOND systems, and 
this clearly indicates that MOND systems are more subject 
to radial-orbit instability than Newtonian models with DM 
and identical total gravitational field. 

In general, for both MOND and ENSs we find a sub- 
stantial spread in the values of ^g, supporting the expecta- 
tion that a "universal" value of does not exist for these 
systems (for Newtonian models with DM, extended versions 
of the stability criterion have been proposed; Polyachenko 
1987; StiaveUi & Sparke 1991). The data in Table 1 suggest 
that in MOND, and even more in the ENSs, increases 
(i.e., relatively more kinetic energy can be stored in radial 
orbits) for decreasing k. While for the ENSs this trend can 
be explained because in the limit k — >■ the stars become 



* Note that by construction a MOND model and its ENS, with 
identical value of ra, have the same value of ^. 



just tracers, it is tempting to speculate that in MOND this 
behaviour can be interpreted as a manifestation of the less 
mixing nature of MOND with respect to Newtonian grav- 
ity. This interpretation is supported by the well established 
numerical finding that phase mixing (Ciotti et al. 2007) and 
violent relaxation (Nipoti et al. 2007a) in MOND systems 
take longer (in units of dynamical times) than in the Newto- 
nian case. Very qualitatively, the deep-MOND force between 
particles behaves like 1/r, i.e. it is nearer to the harmonic 
oscillator force than the force, and it is easy to show 
that a system in which particles interact with the harmonic 
oscillator force, no mixing or instabilities are possible, as 
each particle moves independently of the others (Lynden- 
Bell & Lynden-Bell 1982). Curiously, even though MOND 
forces in a system of particles are non-additive, a similar 
trend towards longer relaxation times has been also found 
recently in shell models interacting with additive 1/r" forces 
(Di Cintio & Ciotti 2011), and the phase-space evolution of 
the case with 1/r forces is strikingly similar to the deep- 
MOND collapses presented in Ciotti et al. (2007). 

The fact that ^ in MOND is strongly affected by the 
properties of the system at large radii suggests to consider 
as a possible stability indicator the quantity ^haif, defined 
as the ratio of radial to tangential kinetic energy within 
''half: the maximum values for stability ^haif.s for our fam- 
ilies of models are reported in the last column of Table 1. 
Remarkably, all our MOND and purely baryonic Newtonian 
families have ^haif.s ~ 1-3, while the ENSs have Chaif,s<;l-5 
(for K = 10) and ^haif,s^3.6 (for k = 0.01), which leads us 
to speculate that a limit on the amount of the radial orbits 
within the half-mass radius might be used as an empirical 
stability criterion for MOND stellar systems, valid from the 
Newtonian to the deep-MOND regime, i.e. independent of 
the value of the internal acceleration relative to oq. This re- 
sult is reminiscent of that of Trenti & Berlin (2006), who 
found that an almost isotropic core can stabilize Newto- 
nian self-gravitating systems with very strong global radial 
anisotropy. 

3.2.3 Simulations of equivalent Newtonian systems with 
live dark-matter halos 

The results on the stability of the ENSs might be affected 
to some extent by our assumption of frozen DM halo, be- 
cause it is possible that an ENS with live DM halo has dif- 
ferent stability properties (see StiaveUi & Sparke 1991). In 
order to assess to what extent our results are affected by 
the assumption of frozen DM halos, we reran the ENS sim- 
ulations of the family EqKio with live DM halos. For tech- 
nical reasons we ran these live-halo simulations with our 
FVFPS treecode (Fortran Version of a Fast Poisson Solver; 
Londrillo, Nipoti, & Ciotti 2003; Nipoti, Londrillo, & Ciotti 
2003), which was already tested against n-mody in Nipoti 
et al. (2007c). As an additional test, we reran with fvfps all 
the purely baryonic Newtonian simulations of the family No, 
finding excellent agreement with those run with n-mody. In 
the live-halo simulations of the models of the family Eqkio 
we used ~ 8 x 10^ particles for the stellar component and 
~ 1.3 X 10® particles for the halo, which has mass ~ 1.6M* 
for the adopted truncation radius rt = lOr* . We verified nu- 
merically that for these two-component EqKio models the 
halo density distribution, which is monotonically decreasing 
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Figure 4. Final axis ratio (c/a)fln as a function of the initial 
anisotropy radius normalized to the half-mass radius for the ENSs 
with 7 = and k = 10 when the DM halo is live (hexagons) and 
frozen (squares). The range of ra/rhaif corresponding to stable 
models (indicated by the arrows) is the same in the two cases. 



(see Fig. 1) and thus satisfies the necessary condition for 
consistency (see Section 2.1), can be generated by a positive 
isotropic distribution function, which we used to distribute 
the DM particles in phase space with the standard rejection 
technique. 

The final axis ratios of the live-halo simulations 
(hexagons in Fig. 4) are not significantly different from 
those of the corresponding frozen-halo simulations (squares 
in Fig. 4). It follows that the stability properties of the Eoftio 
family of models do not depend on whether the halo is mod- 
elled with particles or it is a fixed potential: in particular 
the values of the minimum anisotropy radius ras (and then 
of the maximum Fridman-Polyachenko-Shukhman parame- 
ter ^s) for stability are the same in the two cases. These 
results suggest that our conclusion that MOND models are 
substantially more subject to radial-orbit instability than 
ENSs should be robust, at least if the DM halos are almost 
isotropic. Though the assumption of isotropy of the halos 
is not necessarily justified, given that the anisotropy of DM 
distributions is hard to constrain, it is clear that an isotropic 
halo can be always be invoked, in the context of Newtonian 
gravity, to stabilize an otherwise unstable strongly radially- 
anisotropic stellar system, while this freedom is not allowed 
in the context of MOND. 



4 DISCUSSION AND CONCLUSIONS 

To some extent, the MOND and DM interpretations of the 
kinematics of galaxies can be considered degenerate, i.e. sev- 
eral observational features can be satisfactorily reproduced 
in both paradigms. This raises interesting questions about 
possible tests to discriminate between MOND and Newto- 
nian gravity with DM. Fortunately, it is now well established 



that some important dynamical processes are different, even 
in systems in which the total gravitational potentials are 
identical. Examples are dynamical friction and two-body re- 
laxation. 

In the present work we explored whether radial-orbit in- 
stability acts differently in MOND and in Newtonian grav- 
ity (with and without DM); in addition to the theoretical 
interest, this study can be useful to constrain the interpre- 
tation of the observed kinematics of stellar systems in the 
two cases. In particular, we have focused on the stability 
of OM radially-anisotropic spherical 7-models (7 = and 
7 = 1). We compared the results obtained for MOND mod- 
els, ENSs and purely baryonic Newtonian systems, all with 
the same stellar density distribution. Overall, we found that 
MOND systems are more prone to radial-orbit instability 
than their ENSs, independent of the specific indicator (ra 
or ^) used to quantify the anisotropy. Compared to purely 
baryonic Newtonian systems with the same density profile, 
however, MOND systems have larger minimum anisotropy 
radius for stability r-as, but nevertheless higher maximum 
global anisotropy parameter for stability ^s, a consequence 
of the larger kinetic energy that can be stored in their outer 
regions. We speculate that is larger in MOND systems 
than in Newtonian systems with no DM for the same reasons 
that phase mixing and violent relaxation are less efficient, 
i.e., because the force between particles in deep-MOND de- 
creases with distance less strongly (qualitatively as 1/r) than 
in Newtonian gravity, so being nearer to the special case of 
harmonic oscillator inter-particle forces, when instabilities 
and energy exchanges are impossible 

Observationally, these findings may be relevant to ap- 
plications of MOND to pressure-supported systems such as 
globular clusters, dwarf spheroidal and elliptical galaxies. 
For instance, the velocity-dispersion profiles of globular clus- 
ters in the outer parts of the Milky Way can been used to test 
MOND. In the case of the globular cluster NGC 2419 strong 
OM radial anisotropy (ra close to rac) might be required in 
order to bring the velocity dispersion profile predicted by 
MOND close to the profile inferred from the observed ra- 
dial velocities of the cluster stars (SoUima & Nipoti 2010). 
A combination of the stability constraints discussed in this 
paper with new observations of the radial velocity of stars 
can make the kinematics of NGC 2419 a crucial test for 
MOND (Ibata et al., in preparation). 

Another family of objects that are very interesting from 
the perspective of MOND is that of dwarf spheroidal galax- 
ies, because their low surface-brightness, combined with 
their measured velocity dispersion, lead to conclude that 
they must be DM dominated if Newtonian gravity holds. 
Angus (2008) studied the observed line-of-sight velocity- 
dispersion profiles of Milky Way dwarf spheroidal galaxies, 
finding that in most cases the data can be reproduced in 
MOND, at least with somewhat ad hoc anisotropy profiles. 
This result should be reconsidered on the basis of the con- 
sistency and stability constraints discussed in the present 
work. 

The kinematics of elliptical galaxies has also been con- 
sidered as a possible test for MOND (e.g. Sanders 2000). As 
for globular clusters and dwarf spheroidals, when trying to 
reproduce the kinematics of ellipticals, the anisotropy of the 
velocity distribution of the stars is one of the variables in 
the problem and it is important to constrain it as much as 
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possible, along the lines described in the present paper. For 
instance, consistency constraints (but not the more stringent 
stability constraints) on the anisotropy of stars are consid- 
ered in the recent investigation by Cardone et al. (2010). 
Some authors (Tiret et al. 2007; Klypin & Prada 2009) used 
the observed kinematics of planetary nebulae or of the satel- 
lites around elliptical galaxies as a test for MOND, allowing 
for quite general anisotropy profiles for the satellite system. 
In this case, however, while some limits on the amount of ra- 
dial anisotropy might come from the requirement of consis- 
tency, stability arguments do not necessarily lead to strong 
constraints as long as the considered tracers are not dynam- 
ical dominant. 
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